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Abstract Image distortion due to weak gravitational lensing is examined using a non- 
^s^j ' perturbative method of integrating the geodesic deviation and optical scalar equations 

along the null geodesies connecting the observer to a distant source. The method we 
develop continuously changes the shape of the pencil of rays from the source to the 
observer with no reference to lens planes in astrophysically relevant scenarios. We com- 
pare the projected area and the ratio of semi-major to semi-minor axes of the observed 
elliptical image shape for circular sources from the continuous, thick-lens method with 
the commonly assumed thin-lens approximation. We find that for truncated singular 
isothermal sphere and NFW models of realistic galaxy clusters, the commonly used 
thin-lens approximation is accurate to better than 1 part in 10 4 in predicting the im- 
■ age area and axes ratios. For asymmetric thick lenses consisting of two massive clusters 

separated along the line of sight in redshift up to Az = 0.2, we find that modeling the 
image distortion as two clusters in a single lens plane does not produce relative errors 
in image area or axes ratio more than 0.5%. 
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1 Introduction 



f ^ | The analysis of gravitational lensing systems typically proceeds by invoking a thin- 

lens approximation, where the three-dimensional matter distribution responsible for 
distorting bundles of light rays is projected into a two-dimensional lens plane. Sources 
are assumed to be located in a source plane, and the mapping of image locations (in 
the lens plane) to source locations (in the source plane) is codified into a lens equation 
at the heart of nearly all observational studies. The Jacobian of the lens equation 
map provides the basis for the thin-lens approach to the magnification and shearing of 
images or weak gravitational lensing. 
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Beginning with Frittelli and Newman [TJ and Per lick [2], a number of authors 
developed the theory of gravitational lensing free from lens and source planes by em- 
phasizing that the null geodesic equations provide the exact versions of the "time of 
flight" and lens equations usually invoked in the thin lens approximation. Around the 
same time, Pyne and Birkinshaw _3 approached the problem of light ray paths and 
image magnification by perturbatively solving the null geodesic equations. 

The exact lens equation approach was first applied to spherically symmetric black 
hole solutions [4], [5]. It has since then been applied in the context of lensing by the 
rotating black hole at the center of our galaxy [6], [7]. Other work examines black 
hole lensing, including time delays and magnifications, in the context of relativistic 
images produced by black holes at the centers of galaxies other than our own [§]. For 
cosmological applications widely used in observational strong lensing, the accuracy of 
the thin-lens approximation's predictions of total cluster masses from the appearance 
of arcs or Einstein rings was determined to be approximately 2% in the absence of a 
suitable truncation scheme [5]. 

Following this, a series of papers |10] . [11] redefined image distortion from the exact 
lens equation. The central idea of these papers was to consider the pencil of rays in 
an observer's past light cone that connected the observer to an extended source. The 
cross-sectional shape of the pencil changes along the past light cone. At the observer, 
the projection of the cross-section is the observed shape of the source, while at the 
source, the cross-section of the pencil matches the extended source shape. The shape 
of the cross-section of the pencil was determined by integrating the geodesic deviation 
equations continuously. 

These papers reproduced the formalism of the thin-lens approach by introducing a 
Jacobian of the exact lens mapping. In the first paper [10] . a set of new definitions for 
shape parameters and image distortion parameters were proposed that measured the 
total change of the area or ellipticity (defined as the ratio of semi-major to semi-minor 
axes) of the cross-section of the pencil of rays. These image distortion definitions were 
applied to Schwarzschild black hole geometries. By considering curvature tensors with 
support only in a lens plane, the second paper in the series [TT] reproduced exactly 
the thin-lens Jacobian and discussed the relation between the optical scalars - the 
general relativity "convergence" (or divergence), p, and "shear," a - and the thin- lens 
quantities called the same things: k, and 71 and 72. 

Frittelli and Oberst [T2] then examined image distortion by thick lenses in simple 
(non-physical) models where the Ricci curvature was zero and the Weyl curvature was 
in the form of a square barrier. This paper integrated the optical scalar equations along 
straight line paths through the thick barrier and discussed the distortion parameters 
introduced in [10] . 

From the perspective of the astrophysicist analyzing the statistical, weak-lensing 
signal from real systems, these three papers examining image distortion suffer from 
three significant short-comings. First, the shape and image distortion parameters in- 
troduced in [10] do not correspond in a natural way to the actual observable quantities. 
Second, the curvature models considered in these papers are not appropriate for mod- 
eling galaxy clusters. Finally, the completely general choice of basis vectors (which was 
necessary for examining image distortion in general, say near a black hole) may prevent 
the reader from seeing how one would apply the formalism to typical data where the 
gravitational perturbation is weak. 

Like the work begun in [10] . the purpose of the current paper is to re-examine 
image distortion purely in terms of the physical processes of continuous distortion that 
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occur to a pencil of light rays as the pencil travels from the source to the observer - 
with no reference to lens planes. Unlike the work of Pyne and Birkinshaw, we will solve 
the underlying geodesic and geodesic deviation equations exactly. We will specifically 
assume typical statistical weak lensing conditions, including that the gravitational po- 
tential provides a weak perturbation of a background Robertson- Walker metric. We 
integrate the optical scalar and geodesic deviation equations in this context and state 
our results in the language that is customary in the observational lensing community. 

The standard references, e.g. [13] and [14], provide a very brief and general discus- 
sion of the geodesic deviation equations as the underlying physical description of image 
distortion. This paper explicitly formulates in a direct manner the simplest method of 
integrating the geodesic deviations equations using the Newman-Penrose spin coeffi- 
cient formalism applied in the context of weak gravitational lensing. 

Specifically, we will begin with curvature tensors representing potentially realistic 
matter distributions, and derive expressions for the general relativity convergence, p, 
and shear, a, of pencils of light rays. Then we will integrate the effect of the convergence 
and shear along the pencil from the source to the observer - finding the size and shape 
of the image an observer sees of the object. 

While this approach is generally well known within the general relativity commu- 
nity, its explicit application to image distortion in the case of weak gravitational lensing 
is new. For example, the main equations for image distortion in this method are found 
in [2], but the boundary conditions used there are difficult to apply to the weak lensing 
case because the shape of the pencil is not known a-priori at the observer. 

Because our curvature tensors, convergence, and shear are allowed support through- 
out all of space-time, the approach to image distortion that we will derive represents 
lensing by "thick lenses." The study of thick lenses is becoming increasingly important 
as more examples of merging clusters are found and analyzed, for instance in [15] • The 
physics of merging clusters has important consequences for dark matter models, galaxy 
formation and stellar history, and x-ray modeling. 

Our methods also allow us to begin to examine the accuracy of the thin-lens approx- 
imation in statistical weak gravitational lensing for reasonable lensing configurations. 
We examine the differences in predicted image area and axes ratio for a range of lens 
models. We also perform preliminary investigations of weak gravitational lensing by 
two clusters, closely separated along the line of sight by a varying amount in redshift 
space. 



2 Perturbed metric 

We assume that we have an observer and a series of sources co-moving in a perturbed 
flat Robertson- Walker (RW) space-time. Null geodesies connect the sources to the 
observer and are received simultaneously at the observer. These null geodesies pass 
through a region of space-time where the space-time metric is perturbed by a weak 
gravitational potential, ip. We make the customary assumption that ip depends on 
proper distances at the time the light rays pass the lens. For instance, if the potential 
is spherically symmetric, ip will depend on r p = a(t;)r, where t[ is the time all the null 
geodesies of the observer's past light cone pass the lens. In this paper, as elsewhere, we 
assume a(t;) can be treated as a constant since the scale factor varies slowly compared 
to the time light passes across the lens [13] . The perturbed RW metric is 
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ds 2 = (1 + 2<f) dt 2 - a(t) 2 (l - 2^){dx 2 + dy 2 + dz 2 }, (I) 

where t is the usual co-moving time and (x, y, z) are cartesian coordinates. Throughout 
the paper, we will work in units where G = c = 1. Using the conformal time, the metric, 
Eq. [TJ is conformal to the static metric 

ds 2 = (1 + 2tp) dt 2 - (1 - 2if ){dx 2 + dy 2 + dz 2 }. (2) 

This static metric is a perturbation of a flat, Minkowski-space background metric. 
For pressureless matter, the gravitational potential satisfies 

Vp</5 = AirGpm, (3) 

where p m is a matter density that depends on the proper distance, and the gradient 
operator is taken with respect to the same proper distance: Vp = a(f;) _1 V. 

While the cosmology is "hidden" in Eq. [SJ we will make use of the cosmological 
choice in defining the spatial positions of the lens, source and observer as in [9] . We take 
the background to be a flat RW metric with the Hubble constant Hq — 70 km/s/Mpc, 
the matter density fi m = 0.3 and the cosmological constant density Qj^ = 0.7 = f — Q m - 
In this, the scale factor is |16j 

£2m \ ^ f . ( 3-f/rj yfTTfit^ ^ 



We assume the lens is at a redshift 2; while any sources are at redshifts z s - We use 
the redshift relation, 1 + z = ■ where i e is the time a light ray is emitted that is 

received at an observer at t . Setting a(t ) = 1, we can solve for the value of t e with 
Eq.H 

Then we obtain the radial positions of the source and observer by orienting our 
coordinates in such a way that a light ray travels radially in the background spacetime 
and assuming that the observer, lens and source are at least nearly co- linear. Integrating 
radial null geodesies of the fiat RW metric, Eq. [TJ determines the radial positions of 
the source and observer by ignoring the perturbation introduced by the lens. 



3 Curvature tensors 



In this section, we specify the approximations that apply our formalism to the case of 
ground-based, statistical weak lensing studies. In typical studies, one finds an average 
background redshift of all the sources to compute the projected matter density. This 
is because the measured gravitational shear is related to the dimensionless projected 
matter density k = E/E cr n where £ is the mass density of the lens projected along 
the line of sight and 

E - ° 2Ds 
crlt AitGDiDi, ' 

for angular diameter distances to the lens, source, and between the lens and source. In 
practice, an average background £ C rit, or an average background redshift, is derived 
from averaging color information over all the sources to specify an average D s and Z3; s 
in E cr n. 
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In this paper, we assume that the matter distribution is localized near the origin, 
and that the observer lies at some position along the +z axis. Based on the common 
practice of finding an average background redshift, we take sources to lie at the same 
redshift in a source plane at some fixed value of z < along the —I axis with x 2 + 
y 2 <C z 2 . The observer and sources are far from the region where the gravitational 
perturbation is large. 

The Newman-Penrose (NP) spin coefficient formalism [T7] is based on a tetrad of 
null vectors. Our use of the tetrad vectors is to create complex scalar components of the 
Ricci and Weyl tensors by contracting the curvature tensors with the tetrad vectors. 
The NP tetrad vectors are customarily denoted by 

Aj = {£ ,n ,m ,m ) , (5) 

where the first tetrad vector, £ a , is tangent to the past-directed null geodesic connecting 
the source and observer. The vector m a is a complex spatial vector parallel propagated 
along t a that can be taken to lie in the cross-section of a bundle of light rays surrounding 
t. 

Two NP components of the Ricci and Weyl tensors drive image distortion. The 
Ricci tensor component is 

<%) = —Rabt"^, (6) 
and the Weyl tensor component is 

#0 = -C abcd tm b l c m d . (7) 

In principle, the exact tetrad vectors that are null in the perturbed metric are 
used to find <?oo an d However, we wish to work consistently to first order in the 
gravitational perturbation. Since R a b and C a (,cd are nrs t order in the gravitational 
potential, we must use an approximate tetrad that is null in the flat background metric. 
If (C, C) are complex stereographic coordinates that span a sphere's worth of directions, 
the cartesian coordinate components of the approximate tetrad vectors can be written 



(-l-CC,C + C,*(C-C),-l + CC)> (8) 



\/2(i + CC) 

n a = 1 (l + CC, -(C + C)> i(C - C). 1 - CC) , (9) 

v2(l + CC) 

m ° = V2(i+co (0, 1 " c ~ 2 ' ^ + (10) 

For C = C = 0, the spatial part of t a points parallel to the z axis towards negative 
z values. The approximate tetrad vectors are allowed to vary along the null geodesic 
from the source to the observer, or C(s) where s is a parameter along the geodesic with 
s = at the observer. 

In the context of statistical, wide-field weak gravitational lensing, one studies the 
distortion of source galaxies that lie behind and "near" a large cluster of galaxies that 
serves as a lens. With a typical dithering of the stacked images from a 4-m telescope, 
ground-based studies see approximately 15 arc minutes to each side of the lensing 
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cluster. This means that £ at the observer is small, and, even if it varies along the null 
geodesic, £(s) remains small from the source to the observer. 

For the weakly perturbed metric, Eq.[2j and approximate null tetrad vectors, Eqs.|8] 
and 1101 Kling and Campbell [IS] show that the curvature tensors are related to the 
gravitational potential by 

<%) = ~VV (11) 

and 

3*0 = -X^Pxx ~ Vyy ~ ZiPxy), (12) 

keeping the first order terms of tp to be consistent with the weak-field approximation 
in the metric. All derivatives are taken with respect to the cartesian coordinates, and 
<fixx = d 2 ip/dx 2 . The Ricci tensor, Eq. 111! is exact for all values of (, but the Weyl 
tensor assumes £ to be small and discards terms of order Q(p or smaller. That the Ricci 
and Weyl tensor components do not depend £ implies that these are dependent only on 
the position in space (x(s),y(s), z(s)) along the light ray path for a static perturbation. 



4 Geodesic deviation equations 

The null geodesic connecting the source and observer is part of the observer's past light 
cone. Within this light cone is a pencil of rays that surrounds the null geodesic. The 
pencil of rays collapses to a point at the observer (or the apex of the light cone). The 
shape of this pencil at the source corresponds to the true shape of the source. 

The rate of change of the shape of the pencil is determined by the optical scalars: 
p and a. The divergence, p, measures the rate of the expansion of the pencil. The 
complex shear, a, measures the rate of change of the "shearing" of the pencil, or the 
tendency of the cross section of the pencil to become elliptically shaped. 

If one integrated the null geodesic equations of the perturbed metric, Eq. [2] one 
could find an exact null tetrad {£ a , n a , m a , fh a }. In the spin coefficient formalism, the 
divergence and shear can be computed from the exact null tetrad as 

p = £ a .(,m a m b a = £ a -i ) rn a rn } . (13) 

Alternately, one can compute the optical scalars directly from the curvature tensors 
$00 and 3*o- Using s as the parameter along the null geodesic whose tangent vector is 
£ a , D — -^-^ is the directional derivative along £ a . Then the Sachs equations specify 
how the optical scalars evolve along the null geodesic whose tangent vector is £ a : 

Dp = p 2 + aa + $ 00 , (14) 
Da = 2pcr + !% (15) 

We will treat Eqs. [14] and [15] as ordinary differential equations for p and a where ^oo 
and 3*0 are source terms (functions of s) evaluated at positions (x(s),y(s), z(s)) along 
the null geodesic path. 

Following Penrose and Rindler [T5], we introduce a real Jacobi vector that is carried 
along the null geodesic £ a and is in the cross-section of the pencil of rays: 
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q a = £m a + £m a , (16) 
where m a is the background tetrad vector. Note that since £ is small 

m a ^-L(0,l,-i,0) = -j=(e a x - t e a y ) (17) 

for a pair of constant, real, unit vectors e% and that point in the +x and +y 
directions, respectively. Taking £ = (a — ib)/ \/2, the Jacobi field is related to the real 
basis vectors as 

g a = ae% + fey (18) 
for small Writing Z = (£,£), the geodesic deviation equation is equivalent to 

DZ = -PZ, (19) 

for 



P ° 
a p 



(20) 



Because our bundle of light rays converges at the observer, the divergence, p, is real. 
However, the shear, a is not in general. Writing a = ay + i<Ji, the real and imaginary 
parts of the geodesic deviation equation separate into two real equations: 



Da — — (p + a r )a + (Jib, 

Db = -(p- cr r )b + aid. (21) 

The procedure we will use for integrating the geodesic deviation vectors will be to 
first find the optical scalars from Eqs. [2] and fTS] which are governed by the curvature 
tensors. Then we will integrate the two real components of the geodesic deviation vector 
in Eqs. 1211 The integration of the optical scalars will proceed backwards in time from 
the observer to the source along the past-directed null geodesic from the observer to the 
source. However, because we are ultimately interested in understanding the observed 
shape of the source at the observer, we will integrate the geodesic deviation vectors 
from the source towards the observer (forwards in time). Since we are interested in the 
observational case of the observed shape of intrinsically (on-average) circular sources, 
our approach in integrating the geodesic deviation vectors differs from the presentation 
in typical references in the general relativity community, for example in [2]. 



5 Flat space 

We first consider flat space (ip = 0) to understand the boundary conditions we should 
impose on the integration of p, a and the geodesic deviation vectors. This will also 
motivate our definition of the area and ratio of semi-major to semi-minor axes of 
elliptical images of the source at the observer. 

Here, the curvature tensors are zero, and the spatial part of the null geodesies are 
simply straight lines. Because there is no curvature, the "lens" has axial symmetry 
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about the I axis and we can take a null geodesic whose spatial part is in the x-z plane 
as generic. 

The solutions to the Sachs equations consistent with the pencil forming part of the 
past light cone of the observer are 

p = --!-, . = 0. (22) 

The meaning of the divergence of p at s = is that the rays forming the pencil have 
converged at the observer. 

The geodesic deviation equation, Eg. 1211 becomes 

d = - &=-, (23) 
s s 

where • = d/ds. We set the boundary condition at the source (s = Sf) to be a(sf) = a, 
b(sf) = B with a 2 + /3 2 = 1. The solution to the geodesic deviation equation is then 

as , Bs ,„.. 
a= — b = —. (24) 

s f s f 

We see that the geodesic deviation vector in flat space vanishes at the observer, as it 
should, because the pencil of rays has converged there. 

To model a circular source, we take a pair of geodesic deviation vectors, q a and 
q a , or (a, b) and (a,b), specified by (a = l,/3 = 0) and (a = 0,j3 = 1). For a null 
geodesic whose spatial projection is in the x-z plane, q a points in the +x direction 
and q a points in the +y direction along the entire geodesic, and this pair of geodesic 
deviation vectors uniquely determines the circular cross-section of the bundle along the 
null geodesic from the observer to the source. 

We define the projected ratio of semi-major to semi-minor axes, or axes ratio, to 

be 

R p = lim | -) . (25) 
s^O U I 

and we define the projected area of the image ellipse to be 




In flat space, both the projected area and projected ratio defined this way are 1, which 
means that the image of the circular source that an observer sees is a circle of scaled unit 
area. These definitions for the projected area and ratio hold for any axially symmetric 
lens and will allow for comparison with similarly defined quantities in the traditional 
thin-lens approximation. 



6 Optical scalars with non-zero curvature 



To consider integrating the optical scalar equations with non-zero curvature, we intro- 
duce the variable u = 1/p. The optical scalar equations read 
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ii = -V2(l + u 2 {aa + $ 00 )), (27) 
& = 2V2 (J^J + V2W . (28) 

We assume that the observer is far from any region of significant curvature, so that the 
initial past light cone is a light cone in flat space. The boundary conditions consistent 
with the pencil of rays forming part of the flat-space, past light cone are u(s = 0) = 0, 
a(s — 0) = 0, and 

lim - = 0. 

This last condition insures that if the curvature tensors are zero, one obtains the flat 
space a — solution. These conditions allow for a well defined numerical integration 
of u and a near s — 0. We note from the minus sign in Eq. [23 that u will be negative 
near the observer - so that p ~ — l/(v2s) near the observer. 

At some point after s = 0, it is sometimes convenient to switch back into the (p, a) 
pairing: 



p = V2(p 2 + aa + <%,), (29) 
& = 2V2p<r + \/2>o- (30) 

This is especially true because p can become zero or u — > —00, forcing the switch for 
numerically integrating the equations. 

A general property of pencils of light rays that encounter non-zero curvature at 
some point along the pencil is that conjugate points will form if the geodesic can be 
extended far enough [20]. Conjugate points are places where a Jacobi vector vanishes 
at a value of s 7^ - so that the Jacobi vector vanishes in two places (at the observer 
and somewhere else). 

If the weak energy condition #00 > holds, from Eq. 1291 we see that p is always 
increasing. In fact, p can become zero, meaning that the beam has stopped expanding, 
and then p will eventually run away to positive infinity with a running away to negative 
infinity as the beam contracts. 

Figure [T] shows the behavior in p and a in the case of a point mass curvature for 
two null geodesies: one making an angle of 9 — 60 arc sec from the optical axis and a 
second making an angle of 9 = 120 arc seconds. In the case of the 60 arc sec ray, the 
pencil has reached a conjugate point before s = 0.55. Figure [2] shows the integration 
of two geodesic deviation vector components, (a, 6), backwards in time subject to the 
boundary condition that 

lim — = lim — = 1 

s— >0 S s^0 S 

for s = at the observer. The components are scaled in each case such that a — 1 at 
a final s value. In the case of 9 = 60 arc sec, when the conjugate point is reached, the 
b component of the geodesic deviation vector has vanished. 

The presence of conjugate points presents significant numerical problems to inte- 
grating the optical scalar equations. (The main difficulty is an ambiguity in how to 
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reset the boundary conditions for (p, it) on the other side of the conjugate point.) Con- 
jugate points are common in lensing when one considers rays that pass close to the 
center of a galactic cluster, for instance with Einstein rings. Fortunately for the analysis 
of weak lensing (wide-angles) done on RW cosmological backgrounds with reasonable 
matter distributions representing galaxy clusters, conjugate points tend not to form 
until very, very far beyond the position of an observable lensing source. (The shape of 
a galaxy with redshift more than 1 is very difficult to measure accurately because one 
only detects light from the central region of the galaxy.) 



7 Null geodesic equations 

The optical scalar and geodesic deviation equations, Eqs. 1141 1151 and 1211 are to be 
integrated along the path the light ray takes through the space-time - along the null 
geodesic of the metric, Eq. [2] The simplest way to obtain the null geodesic equations 
is to write out the Euler-Lagrange equations of the Lagrangian 

C- = \g ab i a x = \ ((1 + 2 V ) i 2 -(l- 2<p){x 2 + y 2 + z 2 }) = 0. (31) 

The Lagrangian is zero because the geodesies we seek are null. Working to first order 
in ip and using C = 0, the spatial equations of motion are 

x l = 2x i (1 + 2ip){<p^ k i k ) - 2<f-Vj, (32) 

where tp i = dip/dx 1 . 

For an individual past-directed null geodesies, the boundary conditions at the 
observer, s — 0, are a:(0) = y(0) = 0, z(0) = zq > 0. The initial is deter- 
mined as the radial coordinate distance from the lens (at the origin) when the red- 
shift of the lens is given using Eq. [4] One can set the initial values of (x, y, z) to be 
(sin 9 cos 4>, sin 9 sin <j>, — cos 9) for angles on the observers sky (9,<jj) where 9 = cor- 
responds to looking directly at the lens along the z axis and (j> is a polar coordinate 
around the z axis. (The angle of observation for the center of the image between the 
z axis and the light-ray is 9; the polar angle (j> measures in the "sky plane" an angle 
from the x axis.) 



8 Thin-lens image distortion 

Following Falco, Schnieder and Ehlers [13], the thin-lens approximation involves pro- 
jecting the mass density into a two-dimensional lens plane perpendicular to the optical 
axis and considering a lens mapping from image locations to a source position. The 
lens mapping depends on the bending angle, which is the gradient of the projected 
gravitational potential 

^(x) = — I d 2 x In |x — x'|k(x') 

where the integral is taken over the lens plane. 
The Jacobian of the lens mapping 
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A- 



1 - k - 71 
— 72 1 



-72 
ft + 7i 



(33) 



for the "shears" 71 = (ipxx — ipyy)/2 and 72 = ipxy, controls weak lensing image 
distortion. A^ 1 can be used to map small vectors from the source plane to the lens (or 
image) plane. The eigenvectors of A^ 1 indicate the directions of principle stretching 
and contraction due to lensing, and the eigenvalues give the magnitudes of the lengths 
of the semi-axes. If 7 = \J ~i 2 + 7 2 , then a circular source with unit area is mapped 
into an ellipse where the ratio of the semi-axes is 



Rt., 



1 — K + 7 



1 — K — 7 



(34) 



and the area of the ellipse is 



(1-k) 2 -7 2 



(35) 



9 Axially symmetric examples 

In this section, we compare the thick- lens and thin- lens predictions for image area and 
axes ratio for three axially symmetry matter distributions. We require each matter 
distribution to have a finite total mass when integrated over all space. The matter 
distributions we examine are a point mass model, a singular isothermal sphere (SIS) 
that undergoes hard truncation (matter density set to zero for r > r c ), and a smoothly 
truncated version of the model of Navarro, Frenk, and White (NFW) [21] . For short- 
hand, we will refer to our truncated SIS model as hSIS (reminding the reader of the 
hard truncation) and the smoothly truncated NFW model as tNFW. 

In the next subsection, we briefly outline the three-dimensional potentials and the 
Ricci and Weyl tensor components used in this paper. The formulas for the thin-lens 
projected matter densities, k, and shears, 7, are found in [13], [9], and [21] for the point 
mass, hard truncation SIS, and smoothly truncated NFW. In following subsections, we 
discuss the integration of the optical scalar equations and the appearance of conjugate 
points, and we compare the thin and thick lens image areas and ratios for these axially 
symmetric matter distributions. 



9.1 Thick lens models 



For the point mass model, the gravitational potential is 

* = /2 M 2 2 ' (36) 

where including a; = 1/(1 + 0;) implements the customary approximation that the 

potential depend on the proper distance at the time the light passes the lens [T3] . The 
Ricci tensor, $00 1 is zero an( A the Weyl tensor is 

3M(x - iy) 2 

°~~2a ; (x2+y2 + 2 2)5/2- ( > 
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We take the mass of the point mass model to be 2.0 x 10 solar masses (Mq) 
and place the lens at z/ = 0.45 and sources at z s — 0.8. The cluster mass and redshift 
values chosen here closely correspond to the x-ray luminous cluster RXJ1347-1145, 
and the source redshift corresponds to the average source redshift in one ground-based 
observation of that cluster |22| . 

As discussed in Kling and Frittelli [9], singular isothermal sphere models are un- 
physical without truncation because the total mass of the model (r < oo) is infinite. 
The simplest way to truncate is hard truncation: p — for a proper radius greater than 
some cut-off radius, r c , which is taken as a constant proper distance. The gravitational 
potential is given by 

{1a% In Xp — 2a v 
x p 

where i p = Ojr/r c is a natural dimensionless radial parameter and u v is the velocity 
dispersion - the usual parameter in one-parameter SIS models. At r c , the model has 
been continuously tied to a point mass model with the same total mass. The Ricci 
tensor is 



x p < 1 
x v > 1 



(38) 



for a\r <r c and zero for a;r > r c . The Weyl tensor is 

% = { o 57 ■ ,2 1 1 c < . (40) 

^_ 3r eg „( r ,„) air/rc>1 

We take the cut-off radius to be 3.5 Mpc, which represents a typical size for a 
massive cluster. As in the point mass model, we place the lens at z\ = 0.45, sources at 
z s — 0.8, and set the total mass equal to 2.0 x 1O 15 M0. 

Baltz et al. [21] introduced a smoothly truncated model based on the model of 
Navarro, Frenk, and White (NFW) [23]. The ordinary NFW models, while the best 
fits to n-body simulations of cluster formation at all scales, suffer from unbounded 
total mass in the same way as the un-truncated SIS models. Baltz et al. introduced a 
smooth truncation scheme which makes the total mass bounded by introducing a new 
parameter (the tidal radius, rt) inside of which the NFW and tNFW models agree well. 
We consider the n = 1 truncation scheme, for which matter density is modeled as 

P( r p) = " — — 7T7 ; — ( 41 ) 



m (i+^) 2 (i+Gf) 



Here r p = a;r is a proper radius, p cr ;t is the critical density at the lens redshift, r s is 
a scale radius defined as the peak of r 2 p(r), and 8 C is a characteristic density contrast. 
The characteristic density contrast is related to a concentration parameter, c by 

200 c 3 

Or = 



3 log (1 + c) - c/(l + c) ' 



For our modeling, we take r s = 250 kpc, c = 7.315, and set r = 3c which ensures 
good agreement between the NFW model and the tNFW model within the virial radius 
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[21] . As before, we place the lens at z; = 0.45, sources at z s = 0.8, and the total mass 
for this parameter choice is 2.0 x 10 15 M©. 

The three-dimensional gravitational perturbation is 



GMq 



r s (1 + r 2 ) 5 



tt(t 2 - 1) 



2t 



2 In r + arctan 



{Xp/T)Q; 



T-2- 



-hl 



1 + (xp/r)' 



V 2x p J 



with x = a;r/r s . Here Mo = 47rrs<5 c p cr it and r = rt/r s . <Pqo is given by 



GMpaf r\ 
2r p (r p + r s ) 2 r 2 , + r\ 



The Weyl tensor component is 
9q = 



2{rj + r t ^) 2 



for 



Fi{r P ) = - 



Qr s r1 



2 2 



+ 



4r s 



■raC^pJ = — r + 



rj (r 2 + r 2 ) rf (r 2 + r f ) r$(r$ + rf ) 

2r s 2r s 



2 2 



fp(r p + r s ) rjj(r p + r- s ) 2 r$(r p + r s ) r^(r p + r s ) 2 



(42) 



(43) 



(44) 



3(r t 2 



2r| 



In 



r t( r P + r 's) 2 



Gr s r t 



tan 



9.2 Conjugate points and curvature tensors 

The Ricci and Weyl tensor components #oo an d $b are the source terms in the optical 
scalar equations and are thus the ultimate source to the geodesic deviation equations. 
Understanding the structure of these sources in important to understanding how the 
numerical integration of the geodesic deviation vectors proceeds. 

We have seen that conjugate points can form along the null geodesies, driven by 
divergences in p and a as in Figs. [T] and [2] Whether the optical scalars have divergences 
between the observer and source depends on the size of the perturbation along the given 
geodesic. In general, the unphysical point mass model has higher peaks in the Weyl 
tensor perturbation than more physical models with matter density. 

Figure [3] shows the Ricci and Weyl tensor terms for the three models as a function 
of redshift along a null geodesic in the x-z plane whose angle of observation is 6 — 120 
arc seconds. Both curvature tensors are scaled into dimensionless units. We see that 
compared with the more physical models considered, that the point mass Weyl tensor 
is more strongly peaked, leading to more geodesies with conjugate points. The discrete 
jump in the Ricci tensor term for the hSIS model is due to the matter discontinuity in 
the hard truncation. 
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9.3 Comparisons of thick and thin lens images 

We choose the observer to lie along on the +z axis and the lens to be centered at the 
origin. Because of axial symmetry, choosing the source, lens center, and observer to 
lie in the x-z plane makes the Weyl tensor real along the null geodesic connecting the 
observer and source. Further, a circular source will appear as an ellipse at the observer. 
We restrict ourselves to the typical statistical, weak-lensing observations, where sources 
generally lie outside the critical curve of the lens. For these sources, their projected 
image ellipse at the observer will have the semi-major axis aligned with the y direction. 

For a given model, we can compute the area and ratio of the semi-major to semi- 
minor axis for the observed elliptical image using either the thin-lens version of the 
model and Eqs.[34]and[35] or the ratio and areas from integrating the geodesic deviation 
vectors and using the definitions in Eqs. l25l and l26l We do this by setting an observation 
angle at the observer, 6, and either applying the thin-lens formulas, or tracing a null 
geodesic back in time that makes that initial angle with the z axis at the observer 
until it reaches the source redshift, integrating the optical scalars backwards in time 
along that path. The geodesic deviation equations are integrated forwards in time (from 
source to observer) along the null geodesic using stored values of the optical scalars. 

We can then compute a relative difference between the two predictions: AA/A P — 
1 — A t [/A p and AR/R P = 1 — R t [/R p . Taking the thick-lens values as the true values, 
these two relative differences indicate an error introduced by modeling the system as 
a thin lens. Figures 21 \S\ and [15] show the relative errors in the thin- lens predicted areas 
and axis ratios for the point mass, hSIS, and tNFW models. We see in general that 
the errors are very small. 

In each case, the errors are largest at small radii. The errors generally grow as a 
null geodesies approach the Einstein ring radius, where there is a conjugate point. At 
this conjugate point along the geodesic, and at the caustic in the lens mapping, both 
methods break down. 

The hSIS model considered here has a matter discontinuity at 3.5 Mpc where 
the matter density is (unphysically) set to zero. This corresponds to about 607.8 arc 
sec for our lens and observer positioning. Inside this radius, where there is matter 
density, we see general agreement between the hSIS and tNFW models: the thin-lens 
versions of both models predict smaller axes ratios and areas (positive relative errors) 
for observation angles away from the Einstein ring angle. The dip in the ratio of axes 
in the tNFW model does not appear in the hSIS model. The difference in the two error 
profiles is due to the different internal mass structures of the models. 

Past the truncation radius, the hSIS model becomes the point mass model. Figure[7] 
shows the errors in the ratio and area for the hSIS model near the truncation radius. 
We see that for light rays that do not pass through any of the matter distribution, the 
errors from Fig.[7]are the same as for the point mass model, Fig. [4] small and positive for 
the ratio, small and negative for the area. This result is an internal test of consistency 
for our computational code. The peak in the error for the hSIS model represents rays 
that barely clip the matter distribution. In the case of the smooth, thick-lens approach, 
these rays remain outside the matter distribution except for a very brief segment. Since 
hard truncation is fundamentally unphysical in nature, the peak in the hSIS model is 
not a physical feature, and it should not influence properly designed lensing studies. 
This peak simply indicates an artifact of the unphysical matter modeling. 
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For the most physical model considered, the smoothly truncated NFW, we see that 
the thin-lens approximation very slightly underestimates the area and ratio of elliptical 
images at observation angles away from the Einstein ring angle. 

10 Non-symmetric thick lens 

Next, we consider how to reconstruct elliptical images in the thick lensing approach 
in general. In the axially symmetric case, one knows, a priori, the two specific eigen- 
directions in the cross section of the bundle of light that connects the observer and 
source. For a source, observer, and lens center that are co-planar with the x-z plane, 
these are vectors that point in the +x direction, and +y direction. Note that for a real 
$0, the shear, a, is real, and the geodesic deviation vector equations, Eg. 1211 uncouple. 
If the source lies outside the critical curve in the lens plane, the +x oriented vector will 
shrink in proportion to the +y oriented vector from the source to observer. That the 
geodesic deviation equation uncouples so simply allows one to have a natural, simple, 
basis set for determining shapes. 

In the case of a general lensing configuration, there are still two eigen-directions 
in the bundle of the light from the source to observer that could form a similar basis 
pair to the axially symmetric case. Unfortunately, one does not know these vectors in 
advance. The best method is to follow three geodesic deviation vectors from the source 
to the observer along each null geodesic, and to reconstruct the image ellipse from 
them. 

A general ellipse can be described by 

Y a {t) = Y 1 {t)e% + Y 2 {t)e«, (45) 

with 

Y 1 {t) — L+ cos(i) cos(<5) — L_ sin(t) sin(<5) 

Y 2 (t) = L+ cos(i) sin(<5) + L_ sin(t) cos(<5) (46) 

where L+ and L_ are the semi-major and semi-minor axes of the ellipse, 8 is the angle 
the semi- major axis makes with e%, and t is a parameter that runs from [0, 2tt) |10| . 

If qf — a,ie% + b^ey, = Oj-e£ + bjey, and q^, = a^e% + b^ey are three Jacobi fields, 
denoted collectively by qSj, for each geodesic deviation vector component, we define 

a-, slim- & Ls lim — , (47) 

s^O s s^O s 

to obtain a limiting, projected vector for the image ellipse. Then we have six equations 
for six unknowns that determine the observed, projected ellipse: 

a~y = cos(i-y) cos(<5) — (£— ) sin(fry) sin(5), 

& 7 = (L+) cos(ty) sin(<5) + (L_) sin(ty) cos(<5), (48) 

where ty = (ti,tj,tk) indicate the angle each of the three vectors q® = (qf,qj,q%) 
make along the ellipse measured from the semi-major axis. Therefore, by integrating 
the geodesic deviation equations, Eqs. 1211 simultaneously for three Jacobi fields, one can 
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determine the lengths of the semi- major and semi-minor axes as well as the orientation 
of the ellipse from Eqs. 1481 for general asymmetric lenses. 

As a first test this approach, we model a non-symmetric thick lens by superimposing 
two tNFW profiles whose centers are offset and whose total combined mass is 2 x 
10 15 Mq. We place the observer on the z axis as before and a series of circular sources 
in a source plane at constant — z such that the sources lie at a redshift of z = 0.8. 

The center of the first and slightly more massive cluster, p m \, is positioned on the 
z axis at a redshift of z = 0.448. The center of the second clump, p m 2, lies in the 
x — z plane 20 arc sec from the z axis in the direction of +x. This second clump is at 
a redshift of z = 0.452. With this positioning, the tidal radii of the two clusters barely 
overlap - the center to center proper distance (at a redshift of 0.45) is approximately 
f .1 times the sum of the tidal radii rtx + rt2 (also proper distances). Table Q] indicates 
the parameters used for the model. 

Figure [S] shows what would be the elliptically oriented images of a background sky 
of circular sources. The size of the ellipses drawn here is magnified to show the effects 
of the lensing. Since we only consider objects whose line of sight is greater than 40 arc 
sec away from a position +10 arc sec rightward of the origin, it is reasonable that the 
images are roughly elliptical and not arcs. Also shown are disks corresponding to the 
scale radii of both distributions. Disks for the tidal radii would extend past the borders 
of the figure. One sees the recovery of elliptical shapes, demonstrating the ability to use 
the thick-lens approach to obtain images of the weak lensing of background galaxies. 

11 One lens plane analysis of thick lenses 

Most wide-field surveys used in gravitational lensing detect lensing clusters in portions 
of the sky with no spectroscopic observations, and photometric redshift approximations 
are utilized. Using photometric redshifts only, in the case considered above, where the 
two lensing clusters overlap only at their outmost edges, but are only separated by 
0.004 in redshift, one could not tell that there are two clusters separated along the line 
of sight. Instead, one would model this situation with one lens plane. 

Our thick lens approach to weak lensing allows us to examine the error introduced 
by modeling two separate clusters as being in one lens plane. For the thin lens approx- 
imation, we assume that the projected lensing potential is the sum of two projected 
tNFW model potentials, tjj — tpi + ip-2 , where both lensing potentials are at a distance 
corresponding to a redshift of 0.45. We use the formulas for the tNFW lensing poten- 
tial the appendix of Baltz et al. [21]. We then compute the Jacobian matrix of the 
lens map, and use the eigenvalues and eigenvectors of the thin-lens mapping inverse 
Jacobian matrix to find the shape of the elliptical source following standard references. 

Against the thin lens in a single lens plane, we take the two three-dimensional 
clusters and separate them along the line of sight and integrate the thick-lens shear and 
convergence and Jacobi vectors to find the shape of the image. The observer and center 
of the more massive cluster are always on the +z axis, and the center of the second, 
less massive, cluster is always placed in the x-z plane at an angular separation of 20 arc 
seconds towards +x from the +z axis. The second cluster is always placed at negative 
z values, such that the center of the two clusters (along the z axis) is at a redshift of 
0.45. The redshifts of the two clusters are set according to zx,2 = 0.45 ± i x 0.0004, 
and the greatest separation of redshift values we consider places the two clusters at 
redshifts of 0.44 and 0.46. Photometric redshifts certainly could tell at this greatest 
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separation in redshift that there are two separate clusters. In the modeling for this 
paper, we use the same (c, r s ) pairing as in Table [1] for the two clusters. 

This allows us to plot the relative errors in ellipse ratio and area for the thin- 
lens compared with the thick lens as a function of redshift separation, Az = z<x — z\. 
Figure [9] shows the relative error in the image area, 1 — A-n/Ap, plotted against redshift 
separation for four image locations at (x,y) arc sec separation from origin: (75,0), 
(50,50), (—50,50), (—75,0). Figure [TO] shows the relative error in image axes ratio, 
1 — Rti/Rp, against redshift separation for the same image locations on the sky. In 
both cases, three geodesic deviation vectors were used to first determine the projected 
ellipse seen by the observer, and then the ratio and projected area of the ellipse was 
measured. 

Figures [TT] and [12] show contour plots of the error in the image ellipse's area and 
axes ratio when the two lenses are placed close together, at redshifts of 0.448 and 0.452, 
or far apart, at redshifts of 0.44 and 0.46. 

We see that the value of the relative error in area and axes ratio is positive. In all 
cases, the value of the error at a given point increases as the two clusters are separated. 
However, even at very large separation in redshift, the overall relative error in the axes 
ratio and image area remains small - less than 1% for the area and less than 0.5% for 
the ratio. 



12 Discussion 

Previous papers on gravitational image distortion from the non-perturbative stand- 
point did succeed in laying the ground-work and theory for weak gravitational lensing 
[2], [10], [11], [12]. However, each of these papers examined non- physical lenses, did not 
integrate along null geodesies, and used somewhat opaque language, familiar to the 
relativity community, but unfamiliar to the practicing astrophysics community. This 
paper extends those findings by working through the details of physical lenses, with 
both Ricci and Weyl curvature, and states results in terms of typical, ground-based, 
statistical weak gravitational lensing studies. 

We provide a practical program for integrating the gravitational weak lensing image 
distortion for thick lenses. We draw together the elements of the Newman-Penrose spin- 
coefficient formalism, known to a subset of relativists, and present them in a simplified 
form applicable to weak metric perturbations of a background metric which is conformal 
to the flat Robertson- Walker metric. This program is 

1. place the observer, lens, and sources using measured redshifts, the functional form 
of the scale factor, Eq. [T] and background null geodesies, 

2. simultaneously integrate the null geodesic equations, Eqs. [32] and the Sachs equa- 
tions for the optical scalars, Eqs. 1141 and 1151 from the observer back in time to the 
source, 

3. integrate a combination of geodesic deviation vectors forwards in time from the 
source to the observer using Eq. 1211 and determine the projected shape of the 
image by dividing the components by the differential parameter, s and taking the 
limit as s — > 0. 

We demonstrate this method in a variety of scenarios, including different metric per- 
turbations and non-axially symmetric thick lenses. 
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The traditional weak lensing analysis of a non-axially symmetric lens involves find- 
ing the eigenvectors and eigenvalues of the inverse of the Jacobian of the thin-lens 
mapping. Modern approaches involve solving partial differential equations for the pro- 
jected gravitational potential [21] or the projected mass density [2S]. Integrating the 
true null geodesies has typically been avoided because it was seen as too computation- 
ally intensive. 

The program for finding image distortion from thick lenses that we outline and 
implement numerically in this paper demonstrates that thick-lens analysis has become 
more reasonable. In fact, the time and computationally intensive tasks of gravita- 
tional lensing analysis are image processing, stacking, and shape measurement. While 
thick-lens analysis involving integrating null geodesies and geodesic deviation vectors 
may take longer than the corresponding thin-lens analysis, it is at least an order of 
magnitude smaller task in computation and time than image processing and shape 
measurement. 

Nevertheless, in the case of lensing systems that are axially symmetric, we find that 
the thin-lens approximation generates image ellipses whose area and axes ratio is nearly 
identical to the corresponding thick-lens analysis. The errors for axially symmetric 
systems found in this paper are small compared to observational uncertainties. 

In practice, one observes the averaged ellipticity of sources in a region of the sky 
since no individual background galaxy is known to be circular. For an object with axes 
ratio R, the ellipticity is defined as 

- 1 ~ R 
|£| ~ 1 + R' 

Averaging source ellipticities carries with it an intrinsic error determined by the number 
of objects averaged over and the intrinsic ellipticity dispersion, <r £ , as a = a e /\fN. 
Typical values of <r e are about 0.3. The density of observed images is typically about 
20 useable background objects per square arc minute, as in the DLS survey [26], or 
higher by up to a factor of 2 in other current ground-based studies. 

In analyzing axially symmetric lensing configurations, one typically introduces an- 
nular bins for averaging ellipticity measurements. For instance, one might introduce 
fifteen annular bins of equal width with the center of the first bin at 100 arc sec and 
the center of the last bin at 600 arc sec, and then one would average the ellipticity 
of every object in the bin. Within each bin, there would be an expected number of 
background objects based on the number of objects per square arc minute and the 
size of the annular bin. For a bins at 100 arc sec and 250 arc sec, the observational 
accuracies would be approximately a ~ 0.027 and a « 0.017, respectively, for the DLS 
survey. 

We see that the observational uncertainties are several orders of magnitude larger 
than the error introduced by using a thin-lens approximation. Therefore, the thick 
lens analysis is only of practical benefit if one is attempting to model structure with 
parameters along the line of sight. For example, one might attempt to model a galaxy 
cluster as a triaxial ellipsoid rather than a projection into a two dimensional elliptical 
matter distribution. 

We note that non-truncated SIS models are commonly used in weak gravitational 
lensing studies based on the thin-lens approximation. However, truncation is important 
at the scales of wide-field gravitational lensing analysis - even large clusters are not 
expected to contain a significant mass density past two or three Mpc in radius, which 
is within the field of view for almost all ground-based weak lensing observations of 
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clusters. In this paper, we show that the error in the thin-lens predicted ellipse area 
and axes ratio spikes at the truncation radius. This error spike is a needless artifact of 
poor matter modeling. 

Of significance to the lensing community, is the finding in this paper that a single 
lens plane model of two clusters, separated by a significant distance along the line 
of sight, produces highly accurate image area and axes ratios. Typical weak-lensing 
surveys are not supported by spectroscopic observations of cluster members, but rely 
instead on photometric redshifts to determine cluster membership. We show in this 
paper that the weak lensing image distortion of two clusters that are significantly 
separated in distance along the line of sight but close together on the sky can be 
modeled by a single lens plane to high accuracy. 

This is a problem in two regards. First, studies that attempt to determine the 
number distribution of clusters of different masses at different redshifts that are based 
on gravitational lensing may under-count the total number of clusters, and over-state 
their mass. Second, in studies of structure formation within clusters, if clusters are 
selected on the basis of weak-gravitational lensing image distortion signals, two well- 
separated clusters might be mistaken for a single merging object. 

These findings argue strongly for spectroscopic follow-up to clusters selected from 
weak-gravitational lensing surveys to make sure that the cluster members lie at one 
redshift. 



13 Conclusions 

We have integrated the optical scalar and geodesic deviation equations along the null 
geodesies of a perturbed space-time representing real clusters in a Robertson- Walker 
background. Our calculations demonstrate the feasibility of analyzing thick gravita- 
tional lenses in a non-perturbative scheme and allow an examination of the underlying 
accuracy of the common-thin lens approximation. 

Our principle finding in this context is that the thin-lens approximation is in a sense 
too good. Namely, we find that if there are two lenses significantly separated along the 
line of sight, a single lens plane model of these two lenses reproduces the shapes and 
orientations of image ellipses with an error less than what one could expect from aver- 
aging over the images of close-by uncorrelated background sources. This implies that 
spectroscopic observations are important supplements to weak lensing observations. 

Our initial examinations of the thin-lens approximation for weak-lensing image 
distortion shows that a poor truncation scheme like hard truncation can introduce an 
error that is as large as the inherent thin-lens approximation error in central regions 
of a properly truncated model. Since wide-field, ground-based studies of weak lensing 
do examine the region where the truncation is expected to take place, more careful 
examinations of the influence truncation on the relative error in area and axes ratio is 
warranted. 



14 Appendix: Numerical methods 

In this appendix, we describe details of the numerical calculations in the paper. 

It is most convenient (for numerical integration) to scale all the variables by the age 
of the universe, which we determine by setting a = I in Eq.|4]and solving for t. This has 
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the effect of making distances from the observer to the source have roughly unit size, 
and the parameter s tends to run from to about 0.75 depending on arrangements of 
source and observer. The effect of the lens on the path of the null geodesic and also the 
effect on the optical scalar equations take place over a very short range in s, so adaptive 
step-size methods with error corrections need to be used. In this paper, we follow a 
personal adaptation of the Cash-Karp implementation of the Runge-Kutta-Fehlberg 
4-5 adaptive step-size method |27| . 

The null geodesic path and optical scalars, p (or preferably u = 1/p) and a, are first 
integrated back in time from the observer to the source plane. An initial observation 
angle is set at the observer, and the numerical integration continues until the source 
plane is reached. We refine the step size to stop at the source plane z coordinate value 
with accuracy to one part in 10 12 . As the path and optical scalars are integrated, the 
values of the optical scalars are passed into an array (s, p, a) at every step along the 
numerical integration. 

Since u = 1/p is well defined at the observer, an adaptive step-size algorithm would 
naturally take large steps near the observer. However, p divergences at the observer, 
and this divergence is responsible for making the geodesic deviation vectors vanish at 
the observer. Therefore, we force the adaptive step-size algorithm to take many small 
steps near s = in the integration of the path and optical scalars away from the 
observer towards the source. For most rays, we store approximately 25,000 to 45,000 
steps in the (s, p, a) array, with the vast majority of these steps either near the observer 
or near the lens. 

We then integrate the geodesic deviation equations back from the source to the 
observer, setting the initial conditions at the source to define a circular cross-section 
of the pencil. For this integration, we again use an Runge-Kutta-Fehlberg 4-5 adaptive 
step-size algorithm. We interpolate values for a(s) and p(s) from the stored arrays by 
fitting a fourth order polynomial to the nearest four stored data values. Both optical 
scalars are smooth functions along the integration, although p diverges in the limit 
s -> 0. 

Working to double precision throughout, the total accumulated error in the nu- 
merical integration of null geodesic and optical scalars (from the observer to source) is 
less than one part in 10 14 . The total accumulated error in the numerical integration 
of the geodesic deviation vector components is estimated from the adaptive step-size 
algorithm as less than one part in 10 12 , and the geodesic deviation vector component 
integration typically proceeds in less than 3000 steps. However, further numerical anal- 
ysis of our interpolation indicates that the interpolation might add an error of up to 
one part in 10 10 . 

If s; is the value of the parameter of s along the geodesic at the lens, we find that 
there is essentially no change in the projected components of the geodesic deviation 
vectors, a/s and b/s, for s < sj/20. Therefore, we take the projected limits by setting 
a floor value of s — s t i ny close to s = and stopping the numerical integration there: 

]i m - ~ a ( Stin y ) 

s >Q S ^tiny 

We note that the choice of s t i ny is somewhat arbitrary, but that the best choice is 
related to the maximum step-size set in integrating u near the observer, as well as the 
value of s where the forced small step size ends. 

For the non-axially symmetric cases, we integrate the null geodesic equations and 
optical scalar equations as before, but then simultaneously integrate three geodesic 
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z 



9 (arc sec) 



c 



r a (kpc) Af tot (10 ib M Q ) A 



0.448 
0.452 



0.0 
20.0 



6.43 
6.76 



250 
180 



1.4 
0.6 



1.1 



Table 1 Parameters for two superimposed truncated NFW profile matter distributions in- 
cluding redshift and angle from the z axis of the center of the distribution used for Fig. [8] The 
model parameters and total mass of the cluster (in units of 10 15 solar masses) are also given. 
A = Ar p /(rt\ + rt2) indicates the ratio of the proper distance between the centers of the two 
distributions and the sum of tidal radii. 



deviation vectors back from the source to the observer. The initial conditions for the 
triplet of geodesic deviation vectors are taken such that if the lens was not present the 
triplet of geodesic deviation vectors would be (1,0), (-1/2,^/2), and (l/3,2x/2/3) 
at the observer. This choice of three unrelated vectors insures that Eqs. [48] are not 
linearly dependent. Solving Eqs. [48] is somewhat difficult numerically because of the 
trigonometric functions. Modifying a standard Newton-Raphson method to insure that 
the angular variables do not make large changes (passing 27r, for example) is critical in 
obtaining a solution. Both Maple and Mathematica can find solutions quickly, especially 
if the initial guess is adequate, but both give angular values for ti and Si in unusual 
multiples of n. 

Throughout this paper, error calculations and figures are produced from CH — h 
code based on appropriate algorithms from Numerical Recipes [27]. All calculations 
were independently coded in Mathematica for internal consistency checks. 
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Fig. 2 Geodesic deviation vector components a (solid) and b (dashed) along null geodesies in 
a space-time perturbed by a point mass potential with a total mass of 2 X 10 15 solar masses 
for 9 = 60 and 120 arc sec. The point mass is placed at a redshift of 0.45. The vanishing of the 
b component along the 8 = 60 arc sec null geodesic indicates that a conjugate point has been 
reached. 
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Fig. 3 The Ricci (^oo 5 left) and Weyl ("Z'Oj right) tensor components along a null geodesic 
in the x — z plane that makes an angle of 120 arc sec from the z for point mass (dotted), 
hard-truncated SIS (dashed) and smoothly truncated NFW (solid) models each with a total 
of 2 X 10 15 solar masses. The curvature tensors are scaled into dimensionless units using 
the age of the universe and are plotted against the redshift. The hard-truncation of the SIS 
model shows in the Ricci tensor plot as a step-function-like drop to zero. 
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Fig. 4 The relative error in the ratio of axes (solid) and area (dashed) for a point-mass 
perturbation where the total mass is 2.0 X 10 15 solar masses. 
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Fig. 5 The relative error in the ratio of axes (solid) and area (dashed) for a SIS perturbation 
with hard truncation where the total mass is 2.0 X 10 15 solar masses. The truncation radius 
is set at 3.5 Mpc and the velocity dispersion of the cluster is 1100 km s . The peak in the 
error profile is not well resolved in this figure and is shown in Fig. [7] This peak is located at 
the truncation radius. 
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Fig. 6 The relative error in the ratio of axes (solid) and area (dashed) for a smoothly truncated 
NFW perturbation where the total mass is 2.0 X 10 15 solar masses. The scale radius is set at 
250 kpc, the concentration parameter is set at c = 7.3, and the model's tidal radius is set at 
rt = 3 X c X r s . 
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Fig. 7 The relative error in the ratio of axes (solid) and area (dashed) for the hard truncated 
SIS perturbation in Fig. [5] The hard truncation radius of 3.5 Mpc corresponds to location of the 
peak in the error profile. Outside this radius, the null geodesies see a point mass perturbation. 
This figure shows the consistency in the error plots of the hard truncated SIS and the point 
mass perturbations. 
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Fig. 8 The images of a sky of circular sources lensed by two thick lenses recovered by the 
thick-lens approach. The two clusters are tNFW models at redshifts of 0.448 and 0.452. The 
closer lens is placed on the optical axis and has a total mass of 1.4 X 10 solar masses, while 
the more distant cluster is placed 20 arc sec towards -\-x and has a total mass of 0.6 X 10 15 solar 
masses. The gray and black disks represent the projection of the two clusters scale diameters 
2 X r s on the sky. Each cluster has a tidal diameter extending past 600 arc sec. 
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Fig. 9 The relative error in the projected area as seen by an observer between thick and 
thin-lens analysis for two thick lenses separated in redshift space by Az. The thin-lens analysis 
assumes both lenses are at the same redshift, z = 0.45, which is the center of the redshift 
separation. The two lenses have the same tNFW parameters as in Table [T] The four error 
curves correspond to four different sky image locations: (75, 0) - solid curve, (50, 50) - short 
dashed curve, (—50, 50) - dotted curve, (—75, 0) - long dashed curve. 
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Fig. 10 The relative error in the projected axes ratio as seen by an observer between thick 
and thin-lens analysis for two thick lenses separated in redshift space by Az, following the 
same set-up as Fig. [9] The four error curves correspond to four different sky image locations: 
(75, 0) - solid curve, (50, 50) - short dashed curve, (—50, 50) - dotted curve, (—75, 0) - long 
dashed curve. 
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Fig. 11 A contour plot of the error in the projected area between the thin-lens and thick- 
lens modeling of two lensing clusters. In the left panel, the clusters are positioned redshifts of 
2 = 0.448 and z = 0.452. In the right panel, the redshifts are z = 0.44 and z = 0.46. The disks 
indicate the projection of the size of the scale diameter 2 X r s . 
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Fig. 12 A contour plot of the error in the projected axes ratio between the thin-lens and thick- 
lens modeling of two lensing clusters. In the left panel, the clusters are positioned redshifts of 
z = 0.448 and z = 0.452. In the right panel, the redshifts are z = 0.44 and z = 0.46. The disks 
indicate the projection of the size of the scale diameter 2 X r s . 



